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ABSTRACT 

We investigate the existence of magnetohydrostatic equilibria for topologically complex magnetic fields. 
The approach employed is to perform ideal numerical relaxation experiments. We use a newly-developed 
Lagrangian relaxation scheme that exactly preserves the magnetic field topology during the relaxatio n. Our 
config urations include both twisted and sheared fields, of which some fall into the category for which [Parked 
(1197^ predicted no force-free equilibrium. The first class of field considered contains no magnetic null points, 
and field lines connect between two perfectly conducting plates. In these cases we observe only resolved 
current layers of finite thickness. In further numerical experiments we confirm that magnetic null points are 
loci of singular currents. 

Subject headings: Sun: corona - Sun: magnetic fields 


1. INTRODUCTION 

Magnetic field relaxation in environments like the solar at¬ 
mosphere and laboratory plasma is a crucial process in under¬ 
standing open problems like solar flares and field stability in 
tokamaks. In such environments the field evolves nearly ide¬ 
ally, i.e. the magnetic flux remains frozen to the plasma. For 
an arbitrary br aided inagnetic field between two perfectly con¬ 
ducting planes [Parked (119721) hypothesized that there can be a 
force-free equilibrium of the same topology only if the field’s 
twist varies uniformly along the large-scale magnetic field. 
He further suggested that in resistive magnetohydrodynamics 
(MHD), where reconnection can occur, the field would then 
undergo a rapid change in topology accompanied by magnetic 
energy dissipation that would provide a significant contribu¬ 
tion to coronal heating (lParkedll983bb . 

In subsequent works this idea has been confirmed and chal¬ 
lenged var i ous ti mes ([Parked [1983b[: [Craig & Sneyd[ [2005b 
[Lo^[ 201Q[ [2013[) . Braided magnetic fields from foot point 
motions were shown to be complex enough that they mus t 
exhibit the proposed topological dissipation ([Parked[1983ah . 
[Low[ ([2010[) later showed that there exist solutions for the re¬ 
laxing magnetic field which permit current sheets. One of 
the first siinulation s testing the conjecture was performed by 
[Mikic et H] ([1989h who found filamentary current structures 
with an exponentially increasing strength. Given the limited 
computing power of that time, they were only able to reach 
very moderate resolutions, which renders it questionable if 
they observed proper sheets. 

Doubts about Parke r’s conjecture came from e.g. 
[van BallegooijerJ ([1985h who suggested that a field gen¬ 
erated by foot point motions is able to adjust to those 
motions and reach a force-free state so long as the velocity 
field is continuous at the boundary. This was supported 
by later numerical simulations, in which a series of foot- 
point displacements were performed, and an exponential 
thinning and intensification of current layers was observed 

- rather than a collap se to sub-grid scale of the current 
([van Ballegooijen[[l988h . It has also been suggested that in 
certain configurations no thin current layers - finite o r infinite 

- need necessarily form. [Craig & Snevd[ ([2005h derived 
solutions for relaxing magnetic fields which do not show 
singularities even with sufficiently braided configurations. 


However. [Pontin & Hornig[ ([2015[) recently demonstrated that 
for any braided magnetic field in which the field line mapping 
exhibits small length scales, thin current layers are an in¬ 
evitable feature of the corresponding fo rce-free equilibrium, 
if it ex ists. Building on earlier work bv [Wilmot-Smith et al J 
([2009h they showed that the ideal relaxation of a class of 
braided fields leads to a current distribution of finite strength. 
Moreover, the current layers obtained in the approximate 
force-free equilibria were shown to scale in both thickness 
and intensity with length scales present in t he field line map¬ 
ping, consistent with the earlier results of [van Ballegooijed 

cHi. 

In this work we tackle the problem of current sheet for¬ 
mation during magnetic field relaxation for various topo¬ 
logically non-trivi al configurations at u n preced ented numer¬ 
ical resolution. [Longcope & Straub d 19941) pointed out 
that there exist solutions for relaxed magnetic fields which 
have current layers thin enough that they cannot be distin¬ 
guished from current sheets with moderate grid resolution. 
We apply the newly de veloped numerical code GLEMuR 
([Candelaresi et al.[ [20141) which uses the resources of graph¬ 
ical processing uni ts (GPUs) and makes use of mimetic dif¬ 
ferential operators ([Hvman & Shashko~^[1997h . which greatly 
improve the relaxation quality. The scheme is Lagrangian, 
and is constructed i n such a way that it pe rfectly preserves the 
magnetic topology ([Craig & Snevd[[T9^ . 

Emphasis is put on braids which are not reducible to uni- 
for m twists along a me a n mag netic field such as those used 
by [Wilmot-Smith et all ([20091) . as well a s fields generated 
throug h footpoint motions such as those bv [Longbottom et all 
([1998h . We further investigate the effect of modifying the 
magnetic field to include magnetic nul l points, and show that 
current singularities form there (as in [Pontin & Craig[ [2005b 
[Craig & Pontin[[MT^ . 

2. MODEL AND METHODS 
2.1. Ideal Evolution 

In order to determine existence and structure of equilibria 
for given magnetic topologies, we require to follow an exactly 
ideal evolution. We employ a method that by its construction 
exactly preserves the magnetic fiux, magnetic field line con¬ 
nectivity, and solenoidal nature of the magnetic field B dur¬ 
ing the relaxation. Specifically, we use the Lagrangian code 
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GLEMuR (ICandelaresi et alJ 120141) which solves the equa¬ 
tions for an ideal evolution of a magnetized non-Newtonian 
fluid without inertia, as well as an extension to this method 
that considers a damped fluid with inertia. These methods 
have computational advantages over those that solve for the 
full dynamics of ideal MHD, leading towards a minimum en¬ 
ergy state whose properties are our main concern (rather than 
the evolution to reach the relaxed state). 

In order to preserve the fleld’s topology we make use of 
a Lagrangian grid method where the grid points move along 
with the fluid. If the initial positions of fluid particles at time 
t = 0 are described by the position vector held X, we denote 
their position at time t by x(X, t) with cc(X, 0) = X. These 
fluid elements (grid points) are evolved according to 

( 1 ) 

where the velocity u is chosen in such a way to lead towards 
an equilibrium. We employ different methods for choosing u, 
as outlined below. 

Any ideal evolution of the magnetic held B must be con¬ 
sistent with the ideal induction equation 

BB 

— - V X (-U X B) = 0, (2) 

which implies that the magnetic field is frozen into the fluid 
(lBatcheloirill95Ql: IPriest & Forbesll200Qh . i.e. moves together 
with the fluid particles. From the frozen in condition we can 
relate the magnetic field at later time (following a deformation 
of the fluid particle mesh) to the magnetic field at t = 0; 

= (3) 

j = l ^ 

with Bi being the component of the magnetic field and 
A = det{dxi/dXj) (ICraig & SnevdII 198^: ICandelaresi et al.l 
I2Q14I) . Here the fields are functions of their initial positions 
X and time t. In other words, they are functions of the fluid 
particle positions. 

For son ie of the relaxation sini ulations described herein, 
we follow l^ndelaresi et al.l (1201 4l) by appl ying the magneto- 
frictional term (IChodura & Schlut^[l981h for the evolution 
of the fluid 

u = J X B^ (4) 

with the current density J = X x B. This is the evolu¬ 
tion equation for a non-Newtonian fluid without inertia, and 
the evolution terminates when a force-free field (satisfying 
J X ^ = 0) is attained. This approach is well suited for 
studying relaxation problems, as it is sh own to lead to a mono¬ 
tonic decay of t he magnetic energy (ICraig & Snevdl [198^ 
lYangetaUllQ^ . 

However, there are two disadvantages to this approach. 
First, the monotonic energy decay means that during the re¬ 
laxation the system is unable to escape any small local energy 
minima if a lower global energy minimum exists. Second, in a 
magnetic field containing null points, the null point positions 
are fixed (since the J x B force at the nulls themselves must 
be zero). To address the first issue we consider an extension 
of the method that makes use of inertial effects. The fluid’s 
evolution equation is then given by 

( 5 ) 


with the damping coefficient u and density p. 

To address the second issue of stationary magnetic null 
points we employ a pressure force. In some cases described 
below it is beneficial to seek an equilibrium that is not force- 
free, but where the Forentz force is balanced by a pressure 
gradient. For simplicity here we assume that the pressure is 
directly proportional to the fluid density (corresponding to an 
ideal gas under isothermal changes of state). This yields an 
evolution of the fluid mesh 

u = JxB- pVp, ( 6 ) 

with the compressibility parameter f3. The density can be 
expressed in terms of the initial density po as p{x,t) = 
Po/X = p(X,0)/A, and for convenience we will always 
choose po = 1. We can also add the pressure gradient to 
the inertial evolution equation, to give 

du 

— = {JxB-uu-f3Vp)/p. (7) 

Computing spatial derivatives on a moving grid is a sen¬ 
sitive operation. The direct approach used in previous nu¬ 
merical implementations of the magneto-frictional approach 
involves application of the chain rule leadin g to expressions 
involv ing various products of derivatives (ICraig & Snevdl 
119861) . Using such direct derivatives for computing J = 
V X B on highly distorted grids, such as those we expect 
to occur in our numerical experiments, leads to numerical in¬ 
accuracies, most notably the issue that V • J = 0 is not well 
fulfilled, as was noted bv iPontin et al.l (120091) . Our code GLE¬ 
MuR makes use of mimetic numerical operators to compute 
the curl, which have been shown to more accurately repre¬ 
sent the current on such meshes, and have the advantage that 
they preserve the identity V • (V x B) = 0 up to machine 

S precision for some approp r iate mimetic divergence operator 
Hyman & Shashkovll 19971: ICandelaresi et al.ll2014l) . For the 
time stepping we use a Runge-Kutta 6th order in time ap¬ 
proach. 

All three boundary conditions can be chosen to be periodic 
or line-tied. Here line-tied means that the velocity is set to 
zero and the normal component of the magnet i c field is fixed. 
For studying the problem proposed bv iParkeil (119721) we will 
typically use such line-tied boundaries in the 2 ; direction in 
the simulations described below. But occasionally we will 
impose periodic boundaries. 

2.2. Force-Freeness 

Here we describe some diagnostic tools that are used in 
the following sections to analyse the properties of the final 
states of our relaxation simulations. The evolution of the sys¬ 
tem by equation (|4]) is solely determined by the Forentz force 
Fl = j X B. a force-free state implies Fl = 0, which is 
equivalent to V x F = aB, where a is the force-free param¬ 
eter which satisfies Va-F = 0,i.e.ais constant along mag¬ 
netic field lines. During the relaxation simulations, the mag¬ 
netic field evolves into an energetically more favorable state 
with approximately vanishing Forentz force (when (3 = 0). 
Since the Forentz force never vanishes identically in this nu¬ 
merical approximation, the condition V x F = aB is not 
fulfilled exactly either. We can, nevertheless, still express the 
curl of the magnetic field in terms of a component parallel and 
perpendicular to B : 

V X B = XB - ex B 


du 

dt 


{J X B — uu)Ip, 


( 8 ) 
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with the parameter A and vector e, where we choose e such 
that e B = These two parameters are used to determine 
the deviation from the force-free state quantitatively. 

From equation ([5]) we obtain 


X J B 


’ 

(9) 

J xB 

(10) 

^~B^ 


Comparing A and e for each field line we can infer to what de¬ 
gree the field is force-free. For that we need to trace magnetic 
field lines from the bottom of the domain dXz = Zq to the top 
diiz = Zi and integrate A and \e\ along the lines C: 

X{X,Y) = ^J ^dl, (11) 

c 

c 



c 


where we start our field line integration at (X, Y^Zq). The 
ratio e(X, Y)/X{X^ Y) gives the relative deviation from the 
force-free state. Since S • VA = 0 for the force-free state we 
also comp ute the maximum sl ope of A along the field lines in 
analogy to lPontin et al.l (l2009h and ICandelaresi et aT] (I2Q14I) : 

A*(X,y) = max , (14) 

with the value Ai at point i on the field line and the length of 
the field line 1. 

For magnet i c field lines extending between two parallel 
planes iBerged (119861) suggested a relation between magnetic 
helicity density and the winding of the field around itself. Be¬ 
cause the magnetic helicity density is defined via the mag¬ 
netic vector potential A we choose to measure the twist of the 
magnetic field lines by 

c 


J xB 


B^ 


d/, 


( 12 ) 


For a force-free field this expression reduces to uj{X, Y) = 
s m (q(). _ 

Wilmot-Smith et al] (120101) showed that magnetic field lines 


with a high integrated electric current are places of current 
sheet formation and hence reconnection. In our ideal simu¬ 
lations no reconnection can occur, but of course the forma¬ 
tion of localized current concentrations may take place. To 
analyze their occurrence we compute the magnetic field line 
integrated current density 




dl. 


(16) 


3. BRAIDED FIELDS 

Fro m previous numerical experiments dWilmot-Smith et ^ 
we know that topologically complex braids do not 
necessarily form singular current sheets as the field relaxes 



Fig. 1. — Initial magnetic field lines for the configuration. Colors 
denote the held strength which is strongest in the twist regions. 

towards a force-free state. Here we investigate the re- 
laxation behavior of the m ag netic braids discuss ed by e.g. 
IWilmot-Smith et akl (120091) and lYeates et al.l (120101) . To study 
the relaxation of these fields we use the magneto-frictional 
evolution given by Equation ( 0 ]). 

The initial magnetic field w e consider is the one named 
by I Wilmot-Smith et al.l (120091) . which consists of three braid¬ 
ing regions and a homogeneous background magnetic field 
such that Bz > 0 everywhere. Its form is given by 

B^3 ( 0 ) = Boe^ + - Vc)ex + (x - Xc)ey) 

(17) 

with the initial field strength Bq, strength of twist k, ra¬ 
dius and length in ^-direction of the twist region a and 
I respectively and the twist locations (xc,^c,^c)- We 
choose Xc = {1,-1,1,-1,1,-1}, yc = 0, Zc = 

{-20, -12, -4,4,12,20}, o = \/2, ( = 2 and Bq = 1- To fit 
this configuration into a computational domain, the box size is 
chosen to extend 8 units in x and y and 48 units in 2 ; centered 
at the origin. Upper and lower boundaries are chosen either to 
be line-tied or periodic and the grid resolution is 300 in each 
direction. Sample magnetic field lines are shown in Figure [T] 

3.1. Formation of Current Layers 

As the field evolves and tries to minimize the magnetic en- 
er gy it forms co ncentrations of strong currents. According 
to lParkeil (Il972b singular current sheets should form. How¬ 
ever, we do not find any such formation irrespective of the 
grid resolution (Figure [2l upper panel) and all current concen- 
trations are well resolve d which favors Ballegooijen’s result 
(Ivan Ballegooijenl[l985h . This is even true if we choose peri¬ 
odic boundaries in z-direction (Figure [2l lower panel). 

Varying the grid resolution does not significantly change 
the outcome of these simulations. The width of the current 
layers remains the same, as well as the strength of the current. 

3.2. Topological Complexity 

Since the evolution of the magnetic field is ideal it preserves 
its topology and changes in connectivity are forbidden. One 
measure of the field’s topological complexity is the field line 
integrated electric current density J|| (X, Y). We observe an 
approximate conservation for J 11 for our test configuration of 
(Figure [ 3 ]) . From Figure [3] it is readily seen that despite 
the simple structure of Jz in Figure [2] the thin structure of J|| 
demonstrates the high complexity of the field line configura¬ 
tion E^. 


3.3. Force-Freeness 

Whether or not arbitrarily twisted fiux concentrations are 
allowed to evolve into a force-free state is the second aspect 
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Fig. 2.— z-component of the electric current density at z = 3.8 for times 
close to relaxation for with line tied (upper panel) and periodic boundaries 
in z (lower panel). 

of Parker’s conjecture. Here we monitor the evolution of the 
force-free parameter A*, line averaged Lorentz force e and the 
twist cc for all field lines. 

In line wit h previ ous si mulations bvICraig & Snevdl(ll986b . 
iPontin et ^ (l2009h and ICandelaresi et al l (120141) the field 
evolves such that the domain maximum and average of the 
Lorentz force decreases in time (Figure [4]). This decrease 
is, however, not uniform in the field lines. While e is rather 
smooth at the beginning, it develops large gradients and small- 
scale structures as the field relaxes. In those thin loops the 
Lorentz force no longer decreases and prevents the whole sys¬ 
tem from reaching a force-free equilibrium. 

While e measures the strength of the forces along the field 
lines. A* measures the deviation from the force-free state, i.e. 
V X B = aB with S • Va = 0. As expected, the system 
approaches a state close to force-free (Figure O. At the same 
time it develops small-scale features, like e where A* does not 
change significantly. Those features are a charact eristic of this 
highly twisted field which were illustrated by e.g. lYeates et al.l 
(I201QI) . From Figure [5] we can conclude that, although the 
overall system approaches a more force-free state it does so 
only on average while locally being prevented to reach that 
state. 

By u sing color maps of the magnetic field line lYeates et al.l 
(1201 Ql) showed that regions of different field line mappings are 
connected to a non-trivial topology of the field. Similarly we 


Fig. 3.— Line integrated electric current density J\\(X^ Y) as computed 
from equation for with line tied boundaries at initial time t = 0 and 
time at relaxation f = 60. 

observe regions where the sign of the twist cc changes sharply 
(Figure©. Those are exactly the loci where both e and A* de¬ 
velop into thin structures and A* stays approximately constant 
in time. 

The reason that e and A* develop thin structures as the re¬ 
laxation proceeds is not clear. This could be a feature of 
the numerical method employed to perform the relaxation: 
specifically that under certain conditions the scheme acts to 
reduce the J x ^ force on average within the domain at the 
expense of particular locations at which the relaxation is com¬ 
promised. On the other hand, it is possible that this is asso¬ 
ciated with some more fundamental property of the magnetic 
field. In particular, it could be that the topology of the field, 
as manifested through the sign change of the average field 
line twist cc impedes the further evolution of the field into a 
perfectly force-free state. In order to determine whether this 
is the case we require to develop a theory for the evolution of 
these quantities. This is outside the scope of the present study. 

4. CURRENT FORMATION AT NULL-POINTS 

With the present framework we are able to investigate 
the formation of potentially singular current concentrations 
around magnetic null-points where B = 0. As noted pre¬ 
viously there is strong evidence that in response to appro¬ 
priate perturbations singular current concentrations form at 
nulls in the perfectly conducting limit (ISyrovatskiil 119711: 
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Fig. 4.— Average modulus of the Lorentz force along magnetic field lines 
for the configuration with line tied boundaries at f = 0 (upper panel) and 
t = 60 (lower panel). 

Pontin & Craid 120051: iFuentes-Fernandez & Parnelll 1201 2L 
2013l:ICraig & Pontinll2014l) . Here we embed the null point at 
the base of a coronal loop. In particular, we take the first twist 
region of the magnetic field considered in the previous 
section and insert a parasitic polarity flux patch on the lower 
2 ;-boundary, above which is associated a null point within the 
domain, located at (—0.2229, —0.2229, —7.08330). The sep- 
aratrix surfaces of this null point forms a dome geometry that 
encloses the parasitic polarity. The extent of the domain is 
from (—4, —4, —8) to (4,4, 0) (Figure[7]). 

We study the evolution of this configuration using inertial 
terms and velocity damping and replace the evolution of the 
grid positions by using equation (O. Here we set = 3 and 
choose a grid resolution of 192^. This choice of damping 
term ensures that the magnetic field does not overshoot the 
equilibrium and instead creeps towards it. 

As the field evolves it tries to find a relaxed state of re¬ 
duced Lorentz force. On average over the domain this does 
occur. However, in the absence of plasma pressure, near the 
null point the current density increases to such high values 
that also the Lorentz force starts to diverge, after which the 
simulation stops. The loci of these singular current concentra¬ 
tions are at the magnetic nulls, as is highlighted in Figure [8l 
In line with previous works this current concentration forms 
as the s pine and fan of the null point collapse towards one 
another (iPontin & CraigI 120051: IFuentes-Fernandez & Parnell 


Fig. 5.— Maximum gradient of the force-free parameter A along all mag¬ 
netic field lines for the configuration with line tied boundaries t = 0 
(upper panel) and t = 60 (lower panel). 

nnn. To ensure that this is not a numerical artifact one 
can check that in the absence of the perturbation - i.e. setting 
k = Oin Equation (fTTl) . there is no current growth at the nulls. 
It should be noted that varying the parameter u or resorting to 
the magneto frictional approach does not qualitatively change 
this result. 

Adding a pressure term to our calculations the collapse of 
the fiuid at the magnetic nulls is halted before the numerical 
instability sets in. To achieve this we replace equation © 
by® for the evolution of the fiuid and vary the parameter f3 
which represents the relative weight of the pressure gradient 
to the Lorentz force. Even with the pressure gradient present 
we expect singular current concentrations to form since in 
general the Lorentz force associated with the null point col¬ 
lapse is not irrotation al, and therefore cannot be balanced by 


the pressure gradient ( 

Parnell et al.lll997l:lCrai2 & Litvinenkol 

120051: IPontin & Craid 

20051). Indeed, this is what we observe 


in our simulations where we monitor the maximum current 
II max in the domain at the stage of hydrostatic equilibrium 
(Figure O. By decreasing f3 the maximum current increases, 
as the system gets closer to the zero f3 case. Increasing the grid 
resolution we observe a systematic increase of | J|max, sug¬ 
gesting that we are dealing with a physical cur rent singular¬ 
ity sim ilar to simulations for kink instability by lAli & Snevdl 
(1200 Ih . This holds also true for the case where we replace 
the magneto-frictional term by equation ([5]). As noted by 
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Fig. 6.— Average twist of the magnetic field lines u for the configu¬ 
ration at f = 0 (upper panel) and t = 60 (lower panel). 



Fig. 7.— Initial magnetic field for the configuration with the magnetic 
dome containing the magnetic null and the first twist region of the con¬ 
figuration at the lower boundary at 2 ; = —8. The colors denote the magnetic 
field strength. 

ICraig & Litvinenko! (l2005h : IPontin & Cralgl (120051) . the effect 
of the plasma pressure is to weaken the divergent scaling of 
the peak current density with resolution, indicating that for 
large values of /3 a weaker singularity is present. 

5. SHEARED FIELDS 

Past simulations by iLongbottom et al.l (Il998h of sheared 
magnetic fields suggested the occurrence of singular cur¬ 
rent sheets in the absence of magnetic nulls for sufficiently 



Eig. 8.— Final magnetic field for the first twist region of the E^ con¬ 
figuration with a magnetic dome together with isosurfaces of the magnetic 
field (green, half transparent) and current density (red opaque). For the mag¬ 
netic field we choose a level for the isosurface close to 0 to highlight the area 
around the null, while for J we choose a high value. It can be seen that the 
high current concentration lies at the magnetic null. 



n 


Fig. 9.— Maximum current | J|max at hydrostatic equilibrium for dif¬ 
ferent grid resolutions and pressure parameters /3 for the configuration with 
the magnetic null. The increase with resolution suggests the existence of a 
singular current concentration. 

large shear perturbations. Such fields would then not reach 
a smoo t h forc e-free equilibrium supporting the conjecture of 
iParkeri (ITVtI) . As evidence they pointed to an increasing 
maximum current density as they increased the numerical res¬ 
olution and concluded that the increase will continue indefi¬ 
nitely. As maximum resolution they were able to use 65^ grid 
points. 

Here we propose that their maximum resolution was too 
low to make any meaningful conclusions about the forma¬ 
tion of singular current sheets for cases in which the field 
is highly sheared. As remedy we perform simulations with 
high resolutions and monitor the formation of current lay- 
ers. The field configurati ons are identical to the ones used by 
ILongbottom et al.l (119981) . A Cartesian box of size 2 in each 
dimension is filled with a homogeneous magnetic field in z- 
direction. Subsequently the box is distorted in the ^-direction 
according to 

y = yo- Sa sin {2ttSk{xo + Ox)ILx)z, (18) 
after which we apply a distortion in the x-direction: 

a; = xo - S'b sin(27rS'if(2/+ Oy)/iy) 2 :. (19) 

Here xq and yo are the grid coordinates of the undistorted 
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Fig. 10.— Initial magnetic field for the sheared field configuration with 
Sb = ^ together with the distorted grid box. 

Cartesian grid, Sa and Sb the shearing strengths, Sk the 
wave number, and Oy the origins of the coordinate sys¬ 
tem in X and y and Lx and Ly the length of the box in 
X and y. Here we set the size of the undistorted box to 
Lx = Ly = Lz = 2 and center the domain at the origin. We 
choose Sa = Sk = 1 in all the runs and vary Sb between 
0.1 and 1. Note that the distortion in x-direction is performed 
after the one in ^-direction which is why we use y instead 
of ^0 in Equation ([T9b . For the ^-boundaries we apply the 
line tied condition where the normal component of the field is 
fixed and the grid is rigid. The x and y boundaries are peri¬ 
odic. An example initial configuration is shown in Figure fTOl 
for Sb = 1. 

As the field relaxes towards a more force-free state, the 
maximum current in the simulation domain increases, form¬ 
ing a thin layer running up the center of the domain, centered 
on the z-axis - see Figure [TT] After some time, however, the 
growth of the peak current in the domain flattens off and a 
stable spatio-temporal maximum | J|max is attained. Plotting 
this global maximum of the current as a function of the grid 
resolution for large shears we observe an increase with reso¬ 
lution (Figure [TT]), that eventually saturates - further increase 
of the gird resolution does not lead to an increased current, 
indicating that an underlying finite current layer has been re¬ 
solved. This saturation value strongly increases with shearing 
parameter Sb, as the field distortion prod uces strong currents. 
The reason whv lLongbottom et~n] (Il998h drew the premature 
conclusion that singular current sheets for S'b > 0.4 were 
present was simply due to their limited maximum resolution, 
which suggested that above a certain Sb |^|max would grow 
indefinitely with the resolution suggesting the formation of 
singular current sheets. What is clear is that for the grid res¬ 
olutions they considered an unresolved current concentration 
below the grid scale was present. However, with our high 
resolution simulations we are able to resolve the current con¬ 
centrations even for high grid distortions (Figure [TT]). 

6. CONCLUSIONS 

We have introduced a new computational code (ICandelaresil 
1201 5h that performs an exactly ideal relaxation towards an 
equilibrium magnetic field. This was used to study the prop¬ 
erties of equilibria of various magnetic topologies. Several 
implementations of the relaxation procedure were discussed. 


Fig. 11.— Map of the current density for the sheared field numerical ex¬ 
periments for Sb = 1 at ^ = 0 and final time. Although the current con¬ 
centrates in a small location it is still resolved with the 180^ grid points used 
here and compressed grid cells at high current concentrations. 



Fig. 12.— Maximum current density | J|max in the saturated state in de¬ 
pendence of the grid resolution n for various shearing parameters Sb for the 
sheared field configurations. For all Sb there is eventually a flattening off of 
the curves. 

We implemented both a magneto-frictional approach and an 
approach with velocity damping including plasma inertia. In 
both cases a relaxation towards a force-free state or a magne¬ 
tohydrostatic equilibrium with finite pressure were discussed. 
The code uses a Fagrangian grid approach, and in contrast to 
previous implementations employs mimetic derivatives that 
lead to an improved ap proximation of the final equilibrium 
(ICandelaresi et al.ll201^ . 

We have investigated the ideal evolution of topologically 
non-trivial magnetic field configurations and monitored the 
behavior of the electric current density. The emphasis was on 
determining whether singular current sheets might f orm for 
fields which are sufficiently stressed, as suggested by iParkeil 
(Il972h . Contrary to Parker’s hypothesis we do not find singu¬ 
lar current sheets and all current structures remain resolved in 
the absence of magnetic nulls. 

The first type of field considered was a braided field that 
has been previ ously well studied. In sup port of the previ¬ 
ous results fe.g. lWilmot-Smith et al.ll200'^ we find only well 
resolved current structures. However, we have noted that at 
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contact areas between regions with different field line twist, 
the relaxation of the field towards the force-free state is in¬ 
hibited, as measured by various field line integrated quanti¬ 
ties. This suggests that at least using the artificial path to 
equilibrium discussed here, there may be a barrier to reach¬ 
ing the lowest energy state. This will be discussed further 
in a future publication One should note that, as argued by 
iPontin & HornigI (l2015h . for braided fields of this nature that 
exhibit a field line mapping with very small length scales, 
any equilibrium that does exist must exhibit current layers on 
these same small length scales. Thus, while Parker’s hypothe¬ 
sis for spontaneous formation of current singularities may not 
hold for these fields, the proposal that magnetic braiding can 
provide a source of coronal heating is still valid. In particular, 
as the field is continually braided by the turbulent convective 
motions, the length scales of the current layers will eventually 
become sufficiently small that reconnection occurs. 

We also considered sheared magnetic fields that had previ¬ 
ously been implicated in the formation of current singulari¬ 
ties. We demonstrated that with sufficient grid resolution, a 
finite current layer can always be resolved, in contradiction 
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to the results of iLongbottom et aO (Il998h . who were severely 
limited in the grid resolution available to them. 

Lastly, we have considered magnetic fields containing mag¬ 
netic nulls. We showed that in their presence, strong and un¬ 
resolved current structures form at the ir loci. This has been 
previously observed in variou s studie s (iPontin & Craidl2()0l: 
Fuentes-Fernandez & Parnelll 1201 21 120131: ICraig & PontinI 
20141) . In most of these previous studies a simple linear null 
point was considered. Here we considered a coronal loop with 
a null point near the line-tied boundary in a separatrix dome 
configuration - the perturbation to the field was applied far 
form the null point. Nonetheless, the null point still attracted 
an intense current. 
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